library(limma)
library(DEqMS)
prot_robust <- readRDS('../results/prot_robust.rds')
prot_robust <- readRDS('../results/prot_robust.rds')
notch_per_protein <- readRDS('../results/notch_per_protein.rds')
run_limma <- function(obj, notch_count, name, mock=FALSE){
  if(mock){
    # Treat the first 4 tags as duplicates for two conditions
    extract_cols <- 1:4
    spike <- factor(rep(c('a','b'), c(2,2)))
    contrasts <- 'spikeb'
  } else {
    extract_cols <- 1:10
    spike <- factor(rep(c(1,2,6), c(4,3,3)))
    contrasts <- c('spike2', 'spike6')
  }
  
  contrasts2desc <- c('Mock', '2x vs 1x', '6x vs 1x')
  names(contrasts2desc) <- c('spikeb', 'spike2', 'spike6')
  
  results <- vector('list', length(contrasts))
  names(results) <- unname(contrasts2desc[contrasts])
  
  study.design <- model.matrix(~ spike)
  
  dat <- exprs(obj)[,extract_cols] %>% log(base=2)
  
  fit <- lmFit(dat, study.design)
  fit <- eBayes(fit, trend=TRUE)
  plotSA(fit)
  
  for(contrast in contrasts){

    limma.results <- topTable(fit, coef = contrast, n = Inf)
  
    limma.results <- notch_count %>%
      mutate(sample=remove_x(sample)) %>%
      filter(sample %in% colnames(dat)[extract_cols]) %>%
      group_by(Master.Protein.Accessions) %>%
      summarise(mean_fraction_below=mean(fraction_below)) %>%
      merge(limma.results, by.x='Master.Protein.Accessions', by.y='row.names') %>%
      merge(fData(obj)[,'species',drop=FALSE], by.x='Master.Protein.Accessions', by.y='row.names') %>%
      filter(species!='mixed') %>%
      mutate(contains_notch=mean_fraction_below>0,
             binned_ave_expr=Hmisc::cut2(AveExpr, g=8),
             name=name,
             comparison=contrasts2desc[[contrast]])
    
    results[[contrasts2desc[[contrast]]]] <- limma.results
    
    plot_title <- sprintf('%s - %s', name, contrasts2desc[[contrast]]) 
    p <- ggplot(limma.results) +
      ggtitle(plot_title)
    
    p1 <- p + aes(x = P.Value) + 
      geom_histogram(bins = 50, boundary = 0) + 
      theme_camprot() +
      facet_wrap(~species, scales='free')
    
    print(p1)
    print(p1 + aes(logFC))
    
    p2 <- p +
      aes(log10(P.Value), colour=contains_notch) +
      geom_density() +
      theme_camprot() +
      facet_wrap(~species)
    
    print(p2)
    
    p3 <- p +
      aes(mean_fraction_below, logFC) +
      geom_point() +
      theme_camprot() +
      facet_wrap(~species)
    
    print(p3)
  
    p4 <- p + aes(logFC, colour=contains_notch) +
      geom_density() +
      facet_wrap(~species, scales='free') +
      theme_camprot() +
       xlab('Fold change (Log2)')
  
    if(!mock){
     p4 <- p4 + geom_vline(data=expected[expected$comparison=='2 vs 1',],
                           aes(xintercept=log2(expected)), linetype=2)
   } else{
     p4 <- p4 + geom_vline(xintercept=0 ,linetype=2)
   }
    
    print(p4)
    
    p5 <- limma.results %>%
      group_by(species, contains_notch, binned_ave_expr) %>%
      summarise(proportion_sig=mean(as.numeric(adj.P.Val<0.01))) %>%
      ggplot(aes(binned_ave_expr, proportion_sig, colour=contains_notch, group=contains_notch)) +
      geom_line() +
      theme_camprot(base_size=12) +
      facet_wrap(~species) +
      theme(axis.text.x=element_text(angle=45, vjust=1, hjust=1)) +
      xlab('Average protein abundance') +
      ylab('Propotion significant (1% FDR)') +
      ggtitle(plot_title)

   print(p5)
  }
  
  results <- results %>% do.call(what='rbind')

  results
}

all_mock_limma_results <- summarisations %>% names() %>% lapply(function(method){
    
  summarisations[[method]] %>% names() %>% lapply(function(flt){
    
    summarisations[[method]][[flt]] %>% names() %>% lapply(function(agc){
      obj <- summarisations[[method]][[flt]][[agc]]$protein
      notch_count <- notch_per_protein[[flt]][[agc]]
      
      results <- run_limma(obj, notch_count, mock=TRUE,
                           name=sprintf('%s - %s - %s', agc, flt, method))
      
      results
    }) %>% do.call(what='rbind')
  }) %>% do.call(what='rbind')
}) %>% do.call(what='rbind')%>% separate(name, into=c('agc', 'filtering', 'summarisation'), sep=' - ')
`summarise()` ungrouping output (override with `.groups` argument)

`summarise()` ungrouping output (override with `.groups` argument)

`summarise()` ungrouping output (override with `.groups` argument)

`summarise()` ungrouping output (override with `.groups` argument)

`summarise()` ungrouping output (override with `.groups` argument)

`summarise()` ungrouping output (override with `.groups` argument)

`summarise()` ungrouping output (override with `.groups` argument)

`summarise()` ungrouping output (override with `.groups` argument)

`summarise()` ungrouping output (override with `.groups` argument)

`summarise()` ungrouping output (override with `.groups` argument)

`summarise()` ungrouping output (override with `.groups` argument)



all_limma_results <- summarisations %>% names() %>% lapply(function(method){
    
  summarisations[[method]] %>% names() %>% lapply(function(flt){
    
    summarisations[[method]][[flt]] %>% names() %>% lapply(function(agc){
      obj <- summarisations[[method]][[flt]][[agc]]$protein
      notch_count <- notch_per_protein[[flt]][[agc]]
      
      results <- run_limma(obj, notch_count,
                           name=sprintf('%s - %s - %s', agc, flt, method))
      results
    }) %>% do.call(what='rbind')
  }) %>% do.call(what='rbind')
}) %>% do.call(what='rbind')

`summarise()` ungrouping output (override with `.groups` argument)

`summarise()` ungrouping output (override with `.groups` argument)

`summarise()` ungrouping output (override with `.groups` argument)

`summarise()` ungrouping output (override with `.groups` argument)

`summarise()` ungrouping output (override with `.groups` argument)

`summarise()` ungrouping output (override with `.groups` argument)

`summarise()` ungrouping output (override with `.groups` argument)

`summarise()` ungrouping output (override with `.groups` argument)

`summarise()` ungrouping output (override with `.groups` argument)

`summarise()` ungrouping output (override with `.groups` argument)

`summarise()` ungrouping output (override with `.groups` argument)

NA
NA

for(filt in unique(all_limma_results$filtering)){
  for(method in unique(all_limma_results$summarisation)){  
    p <- all_limma_results %>%
      filter(filtering==filt, summarisation==method) %>%
      group_by(species, contains_notch, agc, comparison) %>%
      summarise(proportion_sig=mean(as.numeric(adj.P.Val<0.01))) %>%
      ggplot(aes(contains_notch, proportion_sig,
                 colour=agc,
                 group=interaction(contains_notch, comparison, agc))) +
      geom_point() +
      theme_camprot(base_size=12) +
      theme(axis.text.x=element_text(angle=45, vjust=1, hjust=1)) +
      xlab('Average protein abundance') +
      ylab('Propotion significant (1% FDR)') +
      facet_grid(comparison~species) +
      ggtitle(sprintf('%s - %s', method, filt)) +
      scale_colour_manual(values=get_cat_palette(2), name='') +
      scale_x_discrete(labels=c('0', '>0'), name='PSMs at/below notch') +
      ylim(0,1)
    
    print(p)

    p <- all_limma_results %>%
      filter(filtering==filt, summarisation==method) %>%
      mutate(binned_ave_expr=as.numeric(Hmisc::cut2(AveExpr, g=4))) %>%
      group_by(species, binned_ave_expr, agc, comparison) %>%
      summarise(proportion_sig=mean(as.numeric(adj.P.Val<0.01), na.rm=TRUE)) %>%
      ggplot(aes(binned_ave_expr, proportion_sig, linetype=agc,
                 group=interaction(comparison, agc))) +
      geom_line() +
      theme_camprot(base_size=12) +
      theme(axis.text.x=element_text(angle=45, vjust=1, hjust=1)) +
      xlab('Average protein abundance') +
      ylab('Propotion significant (1% FDR)') +
      facet_grid(comparison~species) +
      ggtitle(sprintf('%s - %s', method, filt)) +
      scale_linetype_discrete(name='') +
      scale_x_continuous(breaks=1:5, name='Abundance quartile')
    
    print(p)
    
    p <- all_limma_results %>%
      filter(filtering==filt, summarisation==method) %>%
      mutate(binned_ave_expr=as.numeric(Hmisc::cut2(AveExpr, g=4))) %>%
      group_by(species, contains_notch, binned_ave_expr, agc, comparison) %>%
      summarise(proportion_sig=mean(as.numeric(adj.P.Val<0.01), na.rm=TRUE)) %>%
      ggplot(aes(binned_ave_expr, proportion_sig,
                 colour=contains_notch, linetype=agc,
                 group=interaction(contains_notch, comparison, agc))) +
      geom_line() +
      theme_camprot(base_size=12) +
      theme(axis.text.x=element_text(angle=45, vjust=1, hjust=1)) +
      xlab('Average protein abundance') +
      ylab('Propotion significant (1% FDR)') +
      facet_grid(comparison~species) +
      ggtitle(sprintf('%s - %s', method, filt)) +
      scale_colour_manual(values=get_cat_palette(2),
                          name='PSMs at/below notch', labels=c('0', '>0')) +
      scale_linetype_discrete(name='') +
      scale_x_continuous(breaks=1:5, name='Abundance quartile')
    
    print(p)
  }
}
`summarise()` regrouping output by 'species', 'contains_notch', 'agc' (override with `.groups` argument)

print(expected)

for(sp in unique(all_limma_results$species)){
  for(cmp in unique(all_limma_results$comparison)){
    p <- all_limma_results %>%
      filter(species==sp, comparison==cmp) %>%
      group_by(summarisation, contains_notch, agc, filtering) %>%
      summarise(median_fc=median(2^logFC)) %>%
      mutate(filtering=factor(filtering, levels=c('unfiltered', 'filtered', 'filtered, inc S/N'))) %>%
      ggplot(aes(contains_notch, median_fc, colour=summarisation, group=summarisation)) +
      geom_line() +
      theme_camprot(base_size=15) +
      facet_grid(filtering~agc) +
      ggtitle(sprintf('%s - %s', sp, cmp))
    
    print(p)
  }
}

all_limma_results %>%
  filter(summarisation=='sum', filtering=='filtered, inc S/N') %>%
  group_by(contains_notch, agc, filtering, species, comparison) %>%
  summarise(median_fc=median(2^logFC)) %>%
  mutate(comparison=gsub('x', '', comparison)) %>%
  ggplot(aes(contains_notch, median_fc,
             colour=interaction(agc, comparison, sep=':  '),
             group=interaction(agc, comparison, sep=':  '))) +
  geom_line() +
  theme_camprot(base_size=15) +
  facet_wrap(~species) +
  geom_hline(data=expected[expected$comparison %in% c('6 vs 1', '2 vs 1'),],
             aes(yintercept=expected), linetype=2)
for(sp in unique(all_limma_results$species)){
  p <- all_limma_results %>%
    filter(summarisation=='sum', filtering=='filtered, inc S/N', species==sp) %>%
    mutate(comparison=gsub('x', '', comparison)) %>%
    ggplot(aes(logFC,
               linetype=contains_notch,
               colour=comparison)) +
    geom_line(stat='density') +
    theme_camprot(base_size=15) +
    facet_wrap(~agc, scales='free') +
    geom_vline(data=(expected %>% filter(comparison %in% c('6 vs 1', '2 vs 1'), species==sp)),
               aes(xintercept=log2(expected), colour=comparison), linetype=2)
  
  if(sp=='H.sapiens'){
    p <- p + coord_cartesian(xlim=c(-.8,.3))
  }
  print(p)
}

all_mock_limma_results %>%
  group_by(species, filtering, summarisation, contains_notch, agc, comparison) %>%
  summarise(proportion_sig=mean(as.numeric(adj.P.Val<0.01), na.rm=TRUE)) %>%
  arrange(desc(proportion_sig))
`summarise()` regrouping output by 'species', 'filtering', 'summarisation', 'contains_notch', 'agc' (override with `.groups` argument)
all_mock_limma_results %>%
  group_by(species, filtering, summarisation, contains_notch, agc, comparison) %>%
  summarise(median_fc=median(2^logFC)) %>%
  arrange(desc(abs(median_fc))) %>% head()
`summarise()` regrouping output by 'species', 'filtering', 'summarisation', 'contains_notch', 'agc' (override with `.groups` argument)

code for DEqMS if required…

feature_count <- prot_robust$`filtered, inc S/N`$`AGC: 5E4`$feature_counts %>%
  filter(sample %in% colnames(dat))

# Get the minimum peptide count per protein
min.pep.count <- feature_count %>% 
  group_by(Master.Protein.Accessions) %>% 
  summarise(Min.pep.count = min(n)) %>% 
  tibble::column_to_rownames("Master.Protein.Accessions")
  
# Add the min peptide count
fit$count = min.pep.count[rownames(fit$coefficients), "Min.pep.count"]

# Run DEqMS
fit.deqms = spectraCounteBayes(fit)

# Run the DEqMS giagnostic plots
VarianceBoxplot(fit.deqms, n = 30, xlab = "PSM count")
VarianceScatterplot(fit.deqms)

# Extract the DEqMS results
DEqMS.results <- outputResult(fit.deqms, coef_col = 2)
LS0tCnRpdGxlOiAnQWdncmVnYXRlIHRvIHByb3RlaW4tbGV2ZWwgYWJ1bmRhbmNlcycKYXV0aG9yOgogIC0gbmFtZTogIlRvbSBTbWl0aCIKICAgIGFmZmlsaWF0aW9uOiAiQ2FtYnJpZGdlIENlbnRyZSBmb3IgUHJvdGVvbWljcyIKZGF0ZTogImByIGZvcm1hdChTeXMudGltZSgpLCAnJWQgJUIsICVZJylgIgphYnN0cmFjdDogfCAKICBIZXJlLCB3ZSBhZ2dyZWdhdGUgdGhlIFBTTS1sZXZlbCBpbnRlbnNpdGllcyBpbnRvIHByb3RlaW4tbGV2ZWwgaW50ZW5zaXRpZXMKb3V0cHV0OgogIHBkZl9kb2N1bWVudDoKICBodG1sX25vdGVib29rOiBkZWZhdWx0Cmdlb21ldHJ5OiBtYXJnaW49MWluCmZvbnRmYW1pbHk6IG1hdGhwYXpvCmZvbnRzaXplOiAxMXB0Ci0tLQoKYGBge3J9CmxpYnJhcnkoY2FtcHJvdFIpCmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KGxpbW1hKQpsaWJyYXJ5KERFcU1TKQpgYGAKCgpgYGB7cn0KcHJvdF9zdW0gPC0gcmVhZFJEUygnLi4vcmVzdWx0cy9wcm90X3N1bS5yZHMnKQpwcm90X3JvYnVzdCA8LSByZWFkUkRTKCcuLi9yZXN1bHRzL3Byb3Rfcm9idXN0LnJkcycpCm5vdGNoX3Blcl9wcm90ZWluIDwtIHJlYWRSRFMoJy4uL3Jlc3VsdHMvbm90Y2hfcGVyX3Byb3RlaW4ucmRzJykKYGBgCgoKYGBge3J9CnJ1bl9saW1tYSA8LSBmdW5jdGlvbihvYmosIG5vdGNoX2NvdW50LCBuYW1lLCBtb2NrPUZBTFNFKXsKICBpZihtb2NrKXsKICAgICMgVHJlYXQgdGhlIGZpcnN0IDQgdGFncyBhcyBkdXBsaWNhdGVzIGZvciB0d28gY29uZGl0aW9ucwogICAgZXh0cmFjdF9jb2xzIDwtIDE6NAogICAgc3Bpa2UgPC0gZmFjdG9yKHJlcChjKCdhJywnYicpLCBjKDIsMikpKQogICAgY29udHJhc3RzIDwtICdzcGlrZWInCiAgfSBlbHNlIHsKICAgIGV4dHJhY3RfY29scyA8LSAxOjEwCiAgICBzcGlrZSA8LSBmYWN0b3IocmVwKGMoMSwyLDYpLCBjKDQsMywzKSkpCiAgICBjb250cmFzdHMgPC0gYygnc3Bpa2UyJywgJ3NwaWtlNicpCiAgfQogIAogIGNvbnRyYXN0czJkZXNjIDwtIGMoJ01vY2snLCAnMnggdnMgMXgnLCAnNnggdnMgMXgnKQogIG5hbWVzKGNvbnRyYXN0czJkZXNjKSA8LSBjKCdzcGlrZWInLCAnc3Bpa2UyJywgJ3NwaWtlNicpCiAgCiAgcmVzdWx0cyA8LSB2ZWN0b3IoJ2xpc3QnLCBsZW5ndGgoY29udHJhc3RzKSkKICBuYW1lcyhyZXN1bHRzKSA8LSB1bm5hbWUoY29udHJhc3RzMmRlc2NbY29udHJhc3RzXSkKICAKICBzdHVkeS5kZXNpZ24gPC0gbW9kZWwubWF0cml4KH4gc3Bpa2UpCiAgCiAgZGF0IDwtIGV4cHJzKG9iailbLGV4dHJhY3RfY29sc10gJT4lIGxvZyhiYXNlPTIpCiAgCiAgZml0IDwtIGxtRml0KGRhdCwgc3R1ZHkuZGVzaWduKQogIGZpdCA8LSBlQmF5ZXMoZml0LCB0cmVuZD1UUlVFKQogIHBsb3RTQShmaXQpCiAgCiAgZm9yKGNvbnRyYXN0IGluIGNvbnRyYXN0cyl7CgogICAgbGltbWEucmVzdWx0cyA8LSB0b3BUYWJsZShmaXQsIGNvZWYgPSBjb250cmFzdCwgbiA9IEluZikKICAKICAgIGxpbW1hLnJlc3VsdHMgPC0gbm90Y2hfY291bnQgJT4lCiAgICAgIG11dGF0ZShzYW1wbGU9cmVtb3ZlX3goc2FtcGxlKSkgJT4lCiAgICAgIGZpbHRlcihzYW1wbGUgJWluJSBjb2xuYW1lcyhkYXQpW2V4dHJhY3RfY29sc10pICU+JQogICAgICBncm91cF9ieShNYXN0ZXIuUHJvdGVpbi5BY2Nlc3Npb25zKSAlPiUKICAgICAgc3VtbWFyaXNlKG1lYW5fZnJhY3Rpb25fYmVsb3c9bWVhbihmcmFjdGlvbl9iZWxvdykpICU+JQogICAgICBtZXJnZShsaW1tYS5yZXN1bHRzLCBieS54PSdNYXN0ZXIuUHJvdGVpbi5BY2Nlc3Npb25zJywgYnkueT0ncm93Lm5hbWVzJykgJT4lCiAgICAgIG1lcmdlKGZEYXRhKG9iailbLCdzcGVjaWVzJyxkcm9wPUZBTFNFXSwgYnkueD0nTWFzdGVyLlByb3RlaW4uQWNjZXNzaW9ucycsIGJ5Lnk9J3Jvdy5uYW1lcycpICU+JQogICAgICBmaWx0ZXIoc3BlY2llcyE9J21peGVkJykgJT4lCiAgICAgIG11dGF0ZShjb250YWluc19ub3RjaD1tZWFuX2ZyYWN0aW9uX2JlbG93PjAsCiAgICAgICAgICAgICBiaW5uZWRfYXZlX2V4cHI9SG1pc2M6OmN1dDIoQXZlRXhwciwgZz04KSwKICAgICAgICAgICAgIG5hbWU9bmFtZSwKICAgICAgICAgICAgIGNvbXBhcmlzb249Y29udHJhc3RzMmRlc2NbW2NvbnRyYXN0XV0pCiAgICAKICAgIHJlc3VsdHNbW2NvbnRyYXN0czJkZXNjW1tjb250cmFzdF1dXV0gPC0gbGltbWEucmVzdWx0cwogICAgCiAgICBwbG90X3RpdGxlIDwtIHNwcmludGYoJyVzIC0gJXMnLCBuYW1lLCBjb250cmFzdHMyZGVzY1tbY29udHJhc3RdXSkgCiAgICBwIDwtIGdncGxvdChsaW1tYS5yZXN1bHRzKSArCiAgICAgIGdndGl0bGUocGxvdF90aXRsZSkKICAgIAogICAgcDEgPC0gcCArIGFlcyh4ID0gUC5WYWx1ZSkgKyAKICAgICAgZ2VvbV9oaXN0b2dyYW0oYmlucyA9IDUwLCBib3VuZGFyeSA9IDApICsgCiAgICAgIHRoZW1lX2NhbXByb3QoKSArCiAgICAgIGZhY2V0X3dyYXAofnNwZWNpZXMsIHNjYWxlcz0nZnJlZScpCiAgICAKICAgIHByaW50KHAxKQogICAgcHJpbnQocDEgKyBhZXMobG9nRkMpKQogICAgCiAgICBwMiA8LSBwICsKICAgICAgYWVzKGxvZzEwKFAuVmFsdWUpLCBjb2xvdXI9Y29udGFpbnNfbm90Y2gpICsKICAgICAgZ2VvbV9kZW5zaXR5KCkgKwogICAgICB0aGVtZV9jYW1wcm90KCkgKwogICAgICBmYWNldF93cmFwKH5zcGVjaWVzKQogICAgCiAgICBwcmludChwMikKICAgIAogICAgcDMgPC0gcCArCiAgICAgIGFlcyhtZWFuX2ZyYWN0aW9uX2JlbG93LCBsb2dGQykgKwogICAgICBnZW9tX3BvaW50KCkgKwogICAgICB0aGVtZV9jYW1wcm90KCkgKwogICAgICBmYWNldF93cmFwKH5zcGVjaWVzKQogICAgCiAgICBwcmludChwMykKICAKICAgIHA0IDwtIHAgKyBhZXMobG9nRkMsIGNvbG91cj1jb250YWluc19ub3RjaCkgKwogICAgICBnZW9tX2RlbnNpdHkoKSArCiAgICAgIGZhY2V0X3dyYXAofnNwZWNpZXMsIHNjYWxlcz0nZnJlZScpICsKICAgICAgdGhlbWVfY2FtcHJvdCgpICsKICAgICAgIHhsYWIoJ0ZvbGQgY2hhbmdlIChMb2cyKScpCiAgCiAgICBpZighbW9jayl7CiAgICAgcDQgPC0gcDQgKyBnZW9tX3ZsaW5lKGRhdGE9ZXhwZWN0ZWRbZXhwZWN0ZWQkY29tcGFyaXNvbj09JzIgdnMgMScsXSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgYWVzKHhpbnRlcmNlcHQ9bG9nMihleHBlY3RlZCkpLCBsaW5ldHlwZT0yKQogICB9IGVsc2V7CiAgICAgcDQgPC0gcDQgKyBnZW9tX3ZsaW5lKHhpbnRlcmNlcHQ9MCAsbGluZXR5cGU9MikKICAgfQogICAgCiAgICBwcmludChwNCkKICAgIAogICAgcDUgPC0gbGltbWEucmVzdWx0cyAlPiUKICAgICAgZ3JvdXBfYnkoc3BlY2llcywgY29udGFpbnNfbm90Y2gsIGJpbm5lZF9hdmVfZXhwcikgJT4lCiAgICAgIHN1bW1hcmlzZShwcm9wb3J0aW9uX3NpZz1tZWFuKGFzLm51bWVyaWMoYWRqLlAuVmFsPDAuMDEpKSkgJT4lCiAgICAgIGdncGxvdChhZXMoYmlubmVkX2F2ZV9leHByLCBwcm9wb3J0aW9uX3NpZywgY29sb3VyPWNvbnRhaW5zX25vdGNoLCBncm91cD1jb250YWluc19ub3RjaCkpICsKICAgICAgZ2VvbV9saW5lKCkgKwogICAgICB0aGVtZV9jYW1wcm90KGJhc2Vfc2l6ZT0xMikgKwogICAgICBmYWNldF93cmFwKH5zcGVjaWVzKSArCiAgICAgIHRoZW1lKGF4aXMudGV4dC54PWVsZW1lbnRfdGV4dChhbmdsZT00NSwgdmp1c3Q9MSwgaGp1c3Q9MSkpICsKICAgICAgeGxhYignQXZlcmFnZSBwcm90ZWluIGFidW5kYW5jZScpICsKICAgICAgeWxhYignUHJvcG90aW9uIHNpZ25pZmljYW50ICgxJSBGRFIpJykgKwogICAgICBnZ3RpdGxlKHBsb3RfdGl0bGUpCgogICBwcmludChwNSkKICB9CiAgCiAgcmVzdWx0cyA8LSByZXN1bHRzICU+JSBkby5jYWxsKHdoYXQ9J3JiaW5kJykKCiAgcmVzdWx0cwp9CgpgYGAKCgpgYGB7cn0Kc3VtbWFyaXNhdGlvbnMgPC0gbGlzdCgnc3VtJz1wcm90X3N1bSwgJ3JvYnVzdCc9cHJvdF9yb2J1c3QpCgphbGxfbW9ja19saW1tYV9yZXN1bHRzIDwtIHN1bW1hcmlzYXRpb25zICU+JSBuYW1lcygpICU+JSBsYXBwbHkoZnVuY3Rpb24obWV0aG9kKXsKICAgIAogIHN1bW1hcmlzYXRpb25zW1ttZXRob2RdXSAlPiUgbmFtZXMoKSAlPiUgbGFwcGx5KGZ1bmN0aW9uKGZsdCl7CiAgICAKICAgIHN1bW1hcmlzYXRpb25zW1ttZXRob2RdXVtbZmx0XV0gJT4lIG5hbWVzKCkgJT4lIGxhcHBseShmdW5jdGlvbihhZ2MpewogICAgICBvYmogPC0gc3VtbWFyaXNhdGlvbnNbW21ldGhvZF1dW1tmbHRdXVtbYWdjXV0kcHJvdGVpbgogICAgICBub3RjaF9jb3VudCA8LSBub3RjaF9wZXJfcHJvdGVpbltbZmx0XV1bW2FnY11dCiAgICAgIAogICAgICByZXN1bHRzIDwtIHJ1bl9saW1tYShvYmosIG5vdGNoX2NvdW50LCBtb2NrPVRSVUUsCiAgICAgICAgICAgICAgICAgICAgICAgICAgIG5hbWU9c3ByaW50ZignJXMgLSAlcyAtICVzJywgYWdjLCBmbHQsIG1ldGhvZCkpCiAgICAgIAogICAgICByZXN1bHRzCiAgICB9KSAlPiUgZG8uY2FsbCh3aGF0PSdyYmluZCcpCiAgfSkgJT4lIGRvLmNhbGwod2hhdD0ncmJpbmQnKQp9KSAlPiUgZG8uY2FsbCh3aGF0PSdyYmluZCcpICU+JSBzZXBhcmF0ZShuYW1lLCBpbnRvPWMoJ2FnYycsICdmaWx0ZXJpbmcnLCAnc3VtbWFyaXNhdGlvbicpLCBzZXA9JyAtICcpCmBgYAoKYGBge3J9CgoKYWxsX2xpbW1hX3Jlc3VsdHMgPC0gc3VtbWFyaXNhdGlvbnMgJT4lIG5hbWVzKCkgJT4lIGxhcHBseShmdW5jdGlvbihtZXRob2QpewogICAgCiAgc3VtbWFyaXNhdGlvbnNbW21ldGhvZF1dICU+JSBuYW1lcygpICU+JSBsYXBwbHkoZnVuY3Rpb24oZmx0KXsKICAgIAogICAgc3VtbWFyaXNhdGlvbnNbW21ldGhvZF1dW1tmbHRdXSAlPiUgbmFtZXMoKSAlPiUgbGFwcGx5KGZ1bmN0aW9uKGFnYyl7CiAgICAgIG9iaiA8LSBzdW1tYXJpc2F0aW9uc1tbbWV0aG9kXV1bW2ZsdF1dW1thZ2NdXSRwcm90ZWluCiAgICAgIG5vdGNoX2NvdW50IDwtIG5vdGNoX3Blcl9wcm90ZWluW1tmbHRdXVtbYWdjXV0KICAgICAgCiAgICAgIHJlc3VsdHMgPC0gcnVuX2xpbW1hKG9iaiwgbm90Y2hfY291bnQsCiAgICAgICAgICAgICAgICAgICAgICAgICAgIG5hbWU9c3ByaW50ZignJXMgLSAlcyAtICVzJywgYWdjLCBmbHQsIG1ldGhvZCkpCiAgICAgIHJlc3VsdHMKICAgIH0pICU+JSBkby5jYWxsKHdoYXQ9J3JiaW5kJykKICB9KSAlPiUgZG8uY2FsbCh3aGF0PSdyYmluZCcpCn0pICU+JSBkby5jYWxsKHdoYXQ9J3JiaW5kJykgJT4lICU+JSBzZXBhcmF0ZShuYW1lLCBpbnRvPWMoJ2FnYycsICdmaWx0ZXJpbmcnLCAnc3VtbWFyaXNhdGlvbicpLCBzZXA9JyAtICcpCiAgICAgIAoKYGBgCmBgYHtyfQoKZm9yKGZpbHQgaW4gdW5pcXVlKGFsbF9saW1tYV9yZXN1bHRzJGZpbHRlcmluZykpewogIGZvcihtZXRob2QgaW4gdW5pcXVlKGFsbF9saW1tYV9yZXN1bHRzJHN1bW1hcmlzYXRpb24pKXsgIAogICAgcCA8LSBhbGxfbGltbWFfcmVzdWx0cyAlPiUKICAgICAgZmlsdGVyKGZpbHRlcmluZz09ZmlsdCwgc3VtbWFyaXNhdGlvbj09bWV0aG9kKSAlPiUKICAgICAgZ3JvdXBfYnkoc3BlY2llcywgY29udGFpbnNfbm90Y2gsIGFnYywgY29tcGFyaXNvbikgJT4lCiAgICAgIHN1bW1hcmlzZShwcm9wb3J0aW9uX3NpZz1tZWFuKGFzLm51bWVyaWMoYWRqLlAuVmFsPDAuMDEpKSkgJT4lCiAgICAgIGdncGxvdChhZXMoY29udGFpbnNfbm90Y2gsIHByb3BvcnRpb25fc2lnLAogICAgICAgICAgICAgICAgIGNvbG91cj1hZ2MsCiAgICAgICAgICAgICAgICAgZ3JvdXA9aW50ZXJhY3Rpb24oY29udGFpbnNfbm90Y2gsIGNvbXBhcmlzb24sIGFnYykpKSArCiAgICAgIGdlb21fcG9pbnQoKSArCiAgICAgIHRoZW1lX2NhbXByb3QoYmFzZV9zaXplPTEyKSArCiAgICAgIHRoZW1lKGF4aXMudGV4dC54PWVsZW1lbnRfdGV4dChhbmdsZT00NSwgdmp1c3Q9MSwgaGp1c3Q9MSkpICsKICAgICAgeGxhYignQXZlcmFnZSBwcm90ZWluIGFidW5kYW5jZScpICsKICAgICAgeWxhYignUHJvcG90aW9uIHNpZ25pZmljYW50ICgxJSBGRFIpJykgKwogICAgICBmYWNldF9ncmlkKGNvbXBhcmlzb25+c3BlY2llcykgKwogICAgICBnZ3RpdGxlKHNwcmludGYoJyVzIC0gJXMnLCBtZXRob2QsIGZpbHQpKSArCiAgICAgIHNjYWxlX2NvbG91cl9tYW51YWwodmFsdWVzPWdldF9jYXRfcGFsZXR0ZSgyKSwgbmFtZT0nJykgKwogICAgICBzY2FsZV94X2Rpc2NyZXRlKGxhYmVscz1jKCcwJywgJz4wJyksIG5hbWU9J1BTTXMgYXQvYmVsb3cgbm90Y2gnKSArCiAgICAgIHlsaW0oMCwxKQogICAgCiAgICBwcmludChwKQoKICAgIHAgPC0gYWxsX2xpbW1hX3Jlc3VsdHMgJT4lCiAgICAgIGZpbHRlcihmaWx0ZXJpbmc9PWZpbHQsIHN1bW1hcmlzYXRpb249PW1ldGhvZCkgJT4lCiAgICAgIG11dGF0ZShiaW5uZWRfYXZlX2V4cHI9YXMubnVtZXJpYyhIbWlzYzo6Y3V0MihBdmVFeHByLCBnPTQpKSkgJT4lCiAgICAgIGdyb3VwX2J5KHNwZWNpZXMsIGJpbm5lZF9hdmVfZXhwciwgYWdjLCBjb21wYXJpc29uKSAlPiUKICAgICAgc3VtbWFyaXNlKHByb3BvcnRpb25fc2lnPW1lYW4oYXMubnVtZXJpYyhhZGouUC5WYWw8MC4wMSksIG5hLnJtPVRSVUUpKSAlPiUKICAgICAgZ2dwbG90KGFlcyhiaW5uZWRfYXZlX2V4cHIsIHByb3BvcnRpb25fc2lnLCBsaW5ldHlwZT1hZ2MsCiAgICAgICAgICAgICAgICAgZ3JvdXA9aW50ZXJhY3Rpb24oY29tcGFyaXNvbiwgYWdjKSkpICsKICAgICAgZ2VvbV9saW5lKCkgKwogICAgICB0aGVtZV9jYW1wcm90KGJhc2Vfc2l6ZT0xMikgKwogICAgICB0aGVtZShheGlzLnRleHQueD1lbGVtZW50X3RleHQoYW5nbGU9NDUsIHZqdXN0PTEsIGhqdXN0PTEpKSArCiAgICAgIHhsYWIoJ0F2ZXJhZ2UgcHJvdGVpbiBhYnVuZGFuY2UnKSArCiAgICAgIHlsYWIoJ1Byb3BvdGlvbiBzaWduaWZpY2FudCAoMSUgRkRSKScpICsKICAgICAgZmFjZXRfZ3JpZChjb21wYXJpc29ufnNwZWNpZXMpICsKICAgICAgZ2d0aXRsZShzcHJpbnRmKCclcyAtICVzJywgbWV0aG9kLCBmaWx0KSkgKwogICAgICBzY2FsZV9saW5ldHlwZV9kaXNjcmV0ZShuYW1lPScnKSArCiAgICAgIHNjYWxlX3hfY29udGludW91cyhicmVha3M9MTo1LCBuYW1lPSdBYnVuZGFuY2UgcXVhcnRpbGUnKQogICAgCiAgICBwcmludChwKQogICAgCiAgICBwIDwtIGFsbF9saW1tYV9yZXN1bHRzICU+JQogICAgICBmaWx0ZXIoZmlsdGVyaW5nPT1maWx0LCBzdW1tYXJpc2F0aW9uPT1tZXRob2QpICU+JQogICAgICBtdXRhdGUoYmlubmVkX2F2ZV9leHByPWFzLm51bWVyaWMoSG1pc2M6OmN1dDIoQXZlRXhwciwgZz00KSkpICU+JQogICAgICBncm91cF9ieShzcGVjaWVzLCBjb250YWluc19ub3RjaCwgYmlubmVkX2F2ZV9leHByLCBhZ2MsIGNvbXBhcmlzb24pICU+JQogICAgICBzdW1tYXJpc2UocHJvcG9ydGlvbl9zaWc9bWVhbihhcy5udW1lcmljKGFkai5QLlZhbDwwLjAxKSwgbmEucm09VFJVRSkpICU+JQogICAgICBnZ3Bsb3QoYWVzKGJpbm5lZF9hdmVfZXhwciwgcHJvcG9ydGlvbl9zaWcsCiAgICAgICAgICAgICAgICAgY29sb3VyPWNvbnRhaW5zX25vdGNoLCBsaW5ldHlwZT1hZ2MsCiAgICAgICAgICAgICAgICAgZ3JvdXA9aW50ZXJhY3Rpb24oY29udGFpbnNfbm90Y2gsIGNvbXBhcmlzb24sIGFnYykpKSArCiAgICAgIGdlb21fbGluZSgpICsKICAgICAgdGhlbWVfY2FtcHJvdChiYXNlX3NpemU9MTIpICsKICAgICAgdGhlbWUoYXhpcy50ZXh0Lng9ZWxlbWVudF90ZXh0KGFuZ2xlPTQ1LCB2anVzdD0xLCBoanVzdD0xKSkgKwogICAgICB4bGFiKCdBdmVyYWdlIHByb3RlaW4gYWJ1bmRhbmNlJykgKwogICAgICB5bGFiKCdQcm9wb3Rpb24gc2lnbmlmaWNhbnQgKDElIEZEUiknKSArCiAgICAgIGZhY2V0X2dyaWQoY29tcGFyaXNvbn5zcGVjaWVzKSArCiAgICAgIGdndGl0bGUoc3ByaW50ZignJXMgLSAlcycsIG1ldGhvZCwgZmlsdCkpICsKICAgICAgc2NhbGVfY29sb3VyX21hbnVhbCh2YWx1ZXM9Z2V0X2NhdF9wYWxldHRlKDIpLAogICAgICAgICAgICAgICAgICAgICAgICAgIG5hbWU9J1BTTXMgYXQvYmVsb3cgbm90Y2gnLCBsYWJlbHM9YygnMCcsICc+MCcpKSArCiAgICAgIHNjYWxlX2xpbmV0eXBlX2Rpc2NyZXRlKG5hbWU9JycpICsKICAgICAgc2NhbGVfeF9jb250aW51b3VzKGJyZWFrcz0xOjUsIG5hbWU9J0FidW5kYW5jZSBxdWFydGlsZScpCiAgICAKICAgIHByaW50KHApCiAgfQp9CmBgYApgYGB7cn0KcHJpbnQoZXhwZWN0ZWQpCgpmb3Ioc3AgaW4gdW5pcXVlKGFsbF9saW1tYV9yZXN1bHRzJHNwZWNpZXMpKXsKICBmb3IoY21wIGluIHVuaXF1ZShhbGxfbGltbWFfcmVzdWx0cyRjb21wYXJpc29uKSl7CiAgICBwIDwtIGFsbF9saW1tYV9yZXN1bHRzICU+JQogICAgICBmaWx0ZXIoc3BlY2llcz09c3AsIGNvbXBhcmlzb249PWNtcCkgJT4lCiAgICAgIGdyb3VwX2J5KHN1bW1hcmlzYXRpb24sIGNvbnRhaW5zX25vdGNoLCBhZ2MsIGZpbHRlcmluZykgJT4lCiAgICAgIHN1bW1hcmlzZShtZWRpYW5fZmM9bWVkaWFuKDJebG9nRkMpKSAlPiUKICAgICAgbXV0YXRlKGZpbHRlcmluZz1mYWN0b3IoZmlsdGVyaW5nLCBsZXZlbHM9YygndW5maWx0ZXJlZCcsICdmaWx0ZXJlZCcsICdmaWx0ZXJlZCwgaW5jIFMvTicpKSkgJT4lCiAgICAgIGdncGxvdChhZXMoY29udGFpbnNfbm90Y2gsIG1lZGlhbl9mYywgY29sb3VyPXN1bW1hcmlzYXRpb24sIGdyb3VwPXN1bW1hcmlzYXRpb24pKSArCiAgICAgIGdlb21fbGluZSgpICsKICAgICAgdGhlbWVfY2FtcHJvdChiYXNlX3NpemU9MTUpICsKICAgICAgZmFjZXRfZ3JpZChmaWx0ZXJpbmd+YWdjKSArCiAgICAgIGdndGl0bGUoc3ByaW50ZignJXMgLSAlcycsIHNwLCBjbXApKQogICAgCiAgICBwcmludChwKQogIH0KfQoKYWxsX2xpbW1hX3Jlc3VsdHMgJT4lCiAgZmlsdGVyKHN1bW1hcmlzYXRpb249PSdzdW0nLCBmaWx0ZXJpbmc9PSdmaWx0ZXJlZCwgaW5jIFMvTicpICU+JQogIGdyb3VwX2J5KGNvbnRhaW5zX25vdGNoLCBhZ2MsIGZpbHRlcmluZywgc3BlY2llcywgY29tcGFyaXNvbikgJT4lCiAgc3VtbWFyaXNlKG1lZGlhbl9mYz1tZWRpYW4oMl5sb2dGQykpICU+JQogIG11dGF0ZShjb21wYXJpc29uPWdzdWIoJ3gnLCAnJywgY29tcGFyaXNvbikpICU+JQogIGdncGxvdChhZXMoY29udGFpbnNfbm90Y2gsIG1lZGlhbl9mYywKICAgICAgICAgICAgIGNvbG91cj1pbnRlcmFjdGlvbihhZ2MsIGNvbXBhcmlzb24sIHNlcD0nOiAgJyksCiAgICAgICAgICAgICBncm91cD1pbnRlcmFjdGlvbihhZ2MsIGNvbXBhcmlzb24sIHNlcD0nOiAgJykpKSArCiAgZ2VvbV9saW5lKCkgKwogIHRoZW1lX2NhbXByb3QoYmFzZV9zaXplPTE1KSArCiAgZmFjZXRfd3JhcCh+c3BlY2llcykgKwogIGdlb21faGxpbmUoZGF0YT1leHBlY3RlZFtleHBlY3RlZCRjb21wYXJpc29uICVpbiUgYygnNiB2cyAxJywgJzIgdnMgMScpLF0sCiAgICAgICAgICAgICBhZXMoeWludGVyY2VwdD1leHBlY3RlZCksIGxpbmV0eXBlPTIpCmBgYAoKYGBge3J9CmZvcihzcCBpbiB1bmlxdWUoYWxsX2xpbW1hX3Jlc3VsdHMkc3BlY2llcykpewogIHAgPC0gYWxsX2xpbW1hX3Jlc3VsdHMgJT4lCiAgICBmaWx0ZXIoc3VtbWFyaXNhdGlvbj09J3N1bScsIGZpbHRlcmluZz09J2ZpbHRlcmVkLCBpbmMgUy9OJywgc3BlY2llcz09c3ApICU+JQogICAgbXV0YXRlKGNvbXBhcmlzb249Z3N1YigneCcsICcnLCBjb21wYXJpc29uKSkgJT4lCiAgICBnZ3Bsb3QoYWVzKGxvZ0ZDLAogICAgICAgICAgICAgICBsaW5ldHlwZT1jb250YWluc19ub3RjaCwKICAgICAgICAgICAgICAgY29sb3VyPWNvbXBhcmlzb24pKSArCiAgICBnZW9tX2xpbmUoc3RhdD0nZGVuc2l0eScpICsKICAgIHRoZW1lX2NhbXByb3QoYmFzZV9zaXplPTE1KSArCiAgICBmYWNldF93cmFwKH5hZ2MsIHNjYWxlcz0nZnJlZScpICsKICAgIGdlb21fdmxpbmUoZGF0YT0oZXhwZWN0ZWQgJT4lIGZpbHRlcihjb21wYXJpc29uICVpbiUgYygnNiB2cyAxJywgJzIgdnMgMScpLCBzcGVjaWVzPT1zcCkpLAogICAgICAgICAgICAgICBhZXMoeGludGVyY2VwdD1sb2cyKGV4cGVjdGVkKSwgY29sb3VyPWNvbXBhcmlzb24pLCBsaW5ldHlwZT0yKQogIAogIGlmKHNwPT0nSC5zYXBpZW5zJyl7CiAgICBwIDwtIHAgKyBjb29yZF9jYXJ0ZXNpYW4oeGxpbT1jKC0uOCwuMykpCiAgfQogIHByaW50KHApCn0KCmBgYAoKCmBgYHtyfQphbGxfbW9ja19saW1tYV9yZXN1bHRzICU+JQogIGdyb3VwX2J5KHNwZWNpZXMsIGZpbHRlcmluZywgc3VtbWFyaXNhdGlvbiwgY29udGFpbnNfbm90Y2gsIGFnYywgY29tcGFyaXNvbikgJT4lCiAgc3VtbWFyaXNlKHByb3BvcnRpb25fc2lnPW1lYW4oYXMubnVtZXJpYyhhZGouUC5WYWw8MC4wMSksIG5hLnJtPVRSVUUpKSAlPiUKICBhcnJhbmdlKGRlc2MocHJvcG9ydGlvbl9zaWcpKQoKCmFsbF9tb2NrX2xpbW1hX3Jlc3VsdHMgJT4lCiAgZ3JvdXBfYnkoc3BlY2llcywgZmlsdGVyaW5nLCBzdW1tYXJpc2F0aW9uLCBjb250YWluc19ub3RjaCwgYWdjLCBjb21wYXJpc29uKSAlPiUKICBzdW1tYXJpc2UobWVkaWFuX2ZjPW1lZGlhbigyXmxvZ0ZDKSkgJT4lCiAgYXJyYW5nZShkZXNjKGFicyhtZWRpYW5fZmMpKSkgJT4lIGhlYWQoKQoKCgpgYGAKCgoKY29kZSBmb3IgREVxTVMgaWYgcmVxdWlyZWQuLi4KYGBge3J9CmZlYXR1cmVfY291bnQgPC0gcHJvdF9yb2J1c3QkYGZpbHRlcmVkLCBpbmMgUy9OYCRgQUdDOiA1RTRgJGZlYXR1cmVfY291bnRzICU+JQogIGZpbHRlcihzYW1wbGUgJWluJSBjb2xuYW1lcyhkYXQpKQoKIyBHZXQgdGhlIG1pbmltdW0gcGVwdGlkZSBjb3VudCBwZXIgcHJvdGVpbgptaW4ucGVwLmNvdW50IDwtIGZlYXR1cmVfY291bnQgJT4lIAogIGdyb3VwX2J5KE1hc3Rlci5Qcm90ZWluLkFjY2Vzc2lvbnMpICU+JSAKICBzdW1tYXJpc2UoTWluLnBlcC5jb3VudCA9IG1pbihuKSkgJT4lIAogIHRpYmJsZTo6Y29sdW1uX3RvX3Jvd25hbWVzKCJNYXN0ZXIuUHJvdGVpbi5BY2Nlc3Npb25zIikKICAKIyBBZGQgdGhlIG1pbiBwZXB0aWRlIGNvdW50CmZpdCRjb3VudCA9IG1pbi5wZXAuY291bnRbcm93bmFtZXMoZml0JGNvZWZmaWNpZW50cyksICJNaW4ucGVwLmNvdW50Il0KCiMgUnVuIERFcU1TCmZpdC5kZXFtcyA9IHNwZWN0cmFDb3VudGVCYXllcyhmaXQpCgojIFJ1biB0aGUgREVxTVMgZ2lhZ25vc3RpYyBwbG90cwpWYXJpYW5jZUJveHBsb3QoZml0LmRlcW1zLCBuID0gMzAsIHhsYWIgPSAiUFNNIGNvdW50IikKVmFyaWFuY2VTY2F0dGVycGxvdChmaXQuZGVxbXMpCgojIEV4dHJhY3QgdGhlIERFcU1TIHJlc3VsdHMKREVxTVMucmVzdWx0cyA8LSBvdXRwdXRSZXN1bHQoZml0LmRlcW1zLCBjb2VmX2NvbCA9IDIpCgpgYGA=